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Abstract 

Current dissipation in thin superconducting wires is numerically evaluated by 
using the string method, within the framework of time-dependent Ginzburg- 
Landau equation with a Langevin noise term. The most probable transition 
pathway between two neighboring current-carrying metastable states, contin- 
uously linking the Langer-Ambegaokar saddle-point state to a state in which 
the order parameter vanishes somewhere, is found numerically. We also give 
a numerically accurate algorithm to evaluate the prefactors for the rate of 
current-reducing transitions. 
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I. INTRODUCTION 



The picture of resistive (current-reducing) phase slips was first discussed by Little [1]. 
Langer and Ambegaokar then used a Ginzburg-Landau free-energy functional to analyti- 
cally obtain the lowest free-energy saddle point between two current-carrying metastable 
states [2]. The time scale of the resistive phase slips have been formulated by McCum- 
ber and Halperin [3]. The theory developed in Refs. [2] and [3] is generally referred to as 
the Langer-Ambegaokar-McCumber-Halperin (LAMH) theory. Recently, new technique has 
been developed for fabricating superconducting nanowires. Since resistive transition region 
broadens with decreasing cross-sectional area of the wire, nanowires therefore become ideal 
samples for a more precise test of the LAMH theory [4]. For this purpose, a quantita- 
tive evaluation of the thermal-activation rate of phase slip events is needed. In particular, 
this is the case since the thermal rate serves as the background for distinguishing quantum 
fluctuations at low temperature, a topic of considerable basic scientific interest [4] 

Recently, the string method [5-8] has been presented for the numerical evaluation of 
thermally activated rare events. This method first locates the most probable transition 
pathway connecting two metastable states in configuration space. This is done by evolving 
strings, i.e., smooth curves with intrinsic parametrization, into the minimal energy path. 
The transition rates can then be computed using an umbrella sampling technique which 
simulates the fluctuations around the most probable path. In this paper we show that the 
string method can be employed as an efficient numerical tool for the study of thermally 
activated phase slips in thin superconducting wires below T c . 

The system is modeled by a one-dimensional (ID) time- dependent Ginzburg-Landau 
equation (TDGLE) with a Langevin noise term. Applying the string method to this par- 
ticular system, we obtain the most probable transition pathway between two neighbor- 
ing current- carrying metastable states. This pathway continuously connects the Langer- 
Ambegaokar saddle-point state [2] to a state in which the order parameter vanishes some- 
where to allow a phase slip of 2n, as first proposed by Little [1]. We also give a numerically 
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accurate algorithm to evaluate the prefactors for the rate of resistive phase slips. 



II. STRING METHOD 

To outline the string method [5], consider a system governed by the overdamped Langevin 
equation 

75i = -V i y(q) + Ci(0, (!) 

where 7 is the frictional coefficient, q denotes the generalized coordinates {g,}, q\ = dqi/dt, 
Vj = d/dqi, and is a white noise satisfying {0(t)(j(t')) = 2 r ykBT5 i j5(t — £'), with fc B 
denoting the Boltzmann constant and T the temperature. Metastable and stable states are 
located in configuration space as the minima of the potential V(q). Assume q^ and q# are 
the two minima of V. In terms of the topography of V(q), the most probable fluctuation 
which can carry the system from q^ to q B (or q B to q^) corresponds to the lowest intervening 
saddle point between these two minima. The minimal energy path (MEP) is defined as a 
smooth curve q*(s) connecting q^ and q# with intrinsic parametrization such as arc length 
s, which satisfies 

(W) x (q*) = 0, (2) 

where (VV) ± is the component of W normal to the path q*(s). This MEP is the most prob- 
able pathway for thermally activated transitions between q^ and q#. To numerically locate 
the MEP in configuration space, a string q(s) (a smooth curve with intrinsic parametrization 
by s) connecting q^ and q# is evolved according to 

q="(V\/) ± (q). (3) 

A re-parametrization is applied once in a while to enforce accurate parametrization by arc 
length. The stationary solution of Eq. (3) satisfies Eq. (2) which defines the MEP. 

Once the MEP is determined, the lowest saddle point is known and the transition rate 
can be computed by evaluating the fluctuations around the MEP [5]. Following Kramers' 
approach and its generalizations [9-11], the transition rate is given by 
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where q s is the saddle point found at the MEP, H(q) denotes the Hessian of V(q), and \ s is 
the negative eigenvalue of H(q s ). (By definition, H(q s ) has one and only one negative eigen- 
value.) The determinant ratio in Eq. (4) are numerically obtained by linear interpolation 
as follows [7]. 

Let F and G be two N x N positive definite matrices and q a column vector in R N 
space. A harmonic potential parametrized by a (0 < a < 1) is constructed as 
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U a (q) = -q T [(l-a)F + aG]q, 



(5) 



with the corresponding partition function given by 



r(a) = J dq 



exp 



--U a (q) 



(2ire) N / 2 (det[(l - a)F + aG]y l/2 . 
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From the expectation value 
d 
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The expectation value Q(a) can be numerically evaluated in the canonical ensemble governed 
by potential U a (q). In practice, the ensemble is generated by solving 



q = -VU a (q) + <Z, 



(10) 



where ((t) is a white noise satisfying (Ci(t)Cj(t')) = 2e5j-,-5(t — t'). 



To apply the above technique to the present problem, it is noted that the Hessian at 
the saddle point, if(q s ), has a negative eigenvalue \ s . Given this \ s and the corresponding 
normalized eigenvector u s , the indefinite if (q s ) has to be modified to give a positive definite 
#(q s ): 

if (q s ) = H (q s ) + (1/ - A s )u s iiJ, (11) 

where z/ is a positive parameter. It follows that det if (q s ) and det if (q s ) are related by 

det if (q s ) = det if (q s ), 

if we remember that the determinant is the product of the eigenvalues. From the MEP 
q*(s) parametrized by the arc length s, the eigenvector u s can be obtained by evaluating 
dq*(s)/ds at the saddle point, followed by a normalization, and \ s is then computed from 
X s = [u s ] T H(q s )u s . The ratio det if (q^/det if (q s ) can be readily computed according to 
Eq. (9) because if (q^) anc ^ Hfas) are both positive definite. The determinant ratio in the 
rate expression (4) is then obtained from 

det if (q^) v det if (q^) 
detff(q s ) ~ Ys det ff (q,) ' 

III. PHASE-SLIP FLUCTUATIONS IN ONE-DIMENSIONAL 

SUPERCONDUCTOR 

A. One-Dimensional Superconductor 

For a superconducting wire below T c , if the transverse dimension d the coherence 
length £(T), then the variations of the order parameter if) over the cross section of the wire 
are energetically prohibited. The wire sample therefore becomes a ID superconductor, with 
ip being a function of a single coordinate x along the wire. The Ginzburg-Landau free-energy 
functional is of the form 



x)\ = a J dx 



||v^)| 2 - ao(Tc 2 T W )| 2 + ^)| 4 



(12) 



where a is the cross-sectional area of the wire, K = h 2 /m* with m* the effective mass of the 
Cooper pair, and a and (3 are both phenomenological parameters. The time evolution of ijj 
is governed by the time-dependent Ginzburg-Landau equation (TDGLE) 
d 1 SFlib] 

1 dt i) = ~a^T + C = K ^ + a ° (Tc " " m ^ + C ' (13) 
where 7 is a viscosity coefficient, and C( x ?^) is a Langevin white noise, with autocorrelation 
functions 

(C(x, t)C(x', t')) = 0; (C(x, t)C(x', f')> = ^-^ksTSix - x')5(t - f). 

This noise generates a random motion of ip and stabilizes the equilibrium distribution, which 
is proportional to exp{—F[ip(x)]/kBT}. 

For the convenience of presentation and computation, we use the dimensionless form 



dx 



\\vm\ 2 -\\m\ 2 +\\m\ A 



(14) 



for the free-energy functional. Here the over bar denotes the dimensionless quantities, ob- 



tained with F scaled by a^a 2 (T c - T) 2 /{3, ip by ^Ja (T c - T)/ (3, and x by the correlation 



length £(T) = y K/ao(T c — T). Correspondingly, the dimensionless TDGLE is of the form 



= — + c = vv + ^ - w v + c, (15) 

in which the time is scaled by r(T) = ^y/a (T c — T), and the dimensionless noise ( satisfies 
the autocorrelation functions 

(((x,t)((x',t>)) = 0; (C(x,I)C@,F)) = -^^L^ 8{x-x')8{t-t'). 

Throughout the remainder of this paper, all physical quantities are given in terms of the 
dimensionless quantities, using <t£o:q(T c — T) 2 / f3 = crt;H 2 (T)/2ii for the energy scale {a^(T c — 
T) 2 /Aj3 = H 2 (T)/8n is the condensation energy density, where H C (T) is the bulk critical 



field) , \Jao(T c — T)/f3 for the ip scale, £ for the length scale, and r for the time scale. Note 
that all the temperature effects are absorbed into these scales. The over bar will be dropped 
in the remainder of the paper. 



B. Current-Carrying Metastable States 



Consider a closed superconducting ring. Periodic boundary condition for i]){x) is imposed 
by ip(—l/2) = ip(l/2) where / is the circumference of the ring. Metastable current-carrying 
states ip n are obtained from the stationary Ginzburg-Landau equation 

-V 2 V>-V + MV = 0, (16) 

as local minima of F: 

Mx) = fne iknX = yfl - ky k " x , k n = 27in/l (17) 

where k n is the wave vector, f n is the amplitude, and n is an integer. The dimensionless 
current density in the ip n state is given by J n = f%k n = (1 — k^)k n . The metastability of ip n 
requires \k n \ < k c = l/y/3. In the presence of thermodynamic fluctuations, the lifetime of 
these metastable states is finite. When the lifetime is made sufficiently short, the decay of 
persistent current becomes observable. 

C. Current-Reducing Phase-Slip Fluctuations 

The decay of persistent current in a superconducting ring may be explained by dividing 
the ^-function space into different subspaces, each labeled by an integer n, defined through 
4>(l/2) — 4>(— 1/2) = 2im according to the periodic boundary condition for ip = |^|e*^. In 
each ^-function subspace there is a current- carrying metastable state ip n , defined as a local 
minimum of F. 

On the one hand, there are many low-energy configurations that are frequently accessed 
by the fluctuating system. Nevertheless, such low-energy fluctuations cause no change of 
the phase difference (/>(l/2) — <p{—l/2) across the whole ring, and therefore the global phase 
coherence persists. On the other hand, there exist thermodynamic fluctuations that lead to 
transitions between different -^-function subspaces. These fluctuations involve large ampli- 
tude fluctuations of if). Heuristically, if ip vanishes somewhere, then 4>{l/2) — 0(— 1/2) may 
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change (slip) by 2n and hence the system moves from one subspace to another (from ip n 
to VWi or ipn-i) [!]• Since larger persistent current means higher free energy, transitions 
among different metastable states tend to reduce the persistent current on average. Large 
amplitude fluctuations usually cost free energies much higher than ksT, therefore phase-slip 
events are rare. Only when the cross section of the ring is very small and the temperature is 
close enough to T c , the decay of persistent currents due to infrequent phase-slip fluctuations 
becomes observable. 



D. Free-Energy Saddle Point 

Based on the Ginzburg-Landau free-energy functional F, Langer and Ambegaokar have 
derived the lowest saddle point ip s (x) between the two neighboring metastable states ip n (x) 
and ip n -i(x) [2]. This ip s state corresponds to the most probable fluctuation which can 
carry the system from ip n to VVi-i (or from ip n -i to ip n ). Analytical expressions for the 
saddle-point state ip 8 and the dimensionless energy barriers AF_ = F[ip 8 ] — F[ip n ] and 
AF + = F[ip s ] — F[ip n -i] are given in Appendix A. In Sec. V, we will show that using the 
string method, ip s and AF± can be numerically obtained from the MEP connecting ip n (x) 
and VVi-i(z)- 



IV. APPLICATION OF THE STRING METHOD: FLUCTUATION TIME SCALE 



The time scale of the thermodynamic phase-slip fluctuations is determined by the 
Langevin equation (13). The transition rates r± for the transitions ip n -i — > ip n and 
VVi — V'n-i can be written as 

aiHliT) 



r± = fi± exp 



-AF 4 



(18) 



2nk B T 

where f2± are the prefactors which fix the fluctuation time scale, and a^H^(T)/2Tc is the en- 
ergy unit. Based upon Kramers' formulation and its generalizations [9-11], McCumber and 
Halperin have derived an analytical expression for the prefactors Q± [3]. To our knowledge, 
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numerical evaluation of Q± has never been reported. Here we outline a numerical scheme 
for the evaluation of tt± so that a complete solution of the LAMH theory may be obtained. 
Results based on this scheme will be presented in Sec. V. 

We adopt a two-component vector representation for the complex ip: 

rj(x) = [771 (a;) m(x)} T = [Reip(x) Imip(x)] T . 

In terms of rj, the dimensionless form of the free-energy functional in Eq. (14) becomes 



F[n(x)] = J 



dx 



(19) 



and the dimensionless TDGLE becomes 



di 1 = ~6^ F[r1 ^ + C = V 2 T7 + - r/ 2 r7 + C, 
in which the noise = [d(x) (2( X )] T satisfies the autocorrelation functions 



(20) 



(Q(x,t)C 3 (x',t')) 



2k B T 



-5 i:j S{x-x')S(t-t'). 



a € a o(T c — T) 2 /f3 

The vector forms for the metastable state i\) n and the saddle-point state ip s can be easily 
obtained. The Hessian of F[rj(x)] is given by 



H[rj] 



(21) 



-V 2 - 1 + 3r]f + r}l 2rnri2 

277x772 -V 2 - 1 + r]f + 3?7 2 

The general expression (4) for the thermal-activation rate can be directly applied to the 
phase-slip fluctuations, with some elaboration for the symmetry properties of the system. 
According to Eq. (4), the dimensionless form of the prefactors Q± in Eq. (18) can be 
formally written as 



2tt 



det K 



n—l 



det H, 



1/2 



lAi 1 



2tt 



det H n 
I det H, 



1/2 



(22) 



where H n _i, H n , and H s are the three Hessians evaluated at ip n -i, ip n , and ip s according to 
H[rj\ in Eq. (21), and is the lowest (negative) eigenvalue of H s . 
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The free energy F is invariant under the gauge transformation ip(x) — > ip(x)e t ^° and 
the translational transformation ip(x) — > ^(x — rc ). (The gauge invariance in terms of ip is 
equivalent to the rotational invariance in terms of 77). As a consequence, H n must have a zero 
eigenvalue A^ (the lowest one, other eigenvalues are all positive), with the corresponding 
eigenvector 



ui 1} (*) 



a 



(i) 







1 



(23) 



where r) n is the vector form for ip n and is the normalization factor. Similarly, H s must 
have a zero eigenvalue \ f\ with the corresponding eigenvector 



uf\x) = a 



(2) 
s 



-1 

1 



(24) 



where r} s is the vector form for ip s and a^f 1 is the normalization factor. The translational 
invariance of F leads to another zero eigenvalue A® because the phase-slip center in the 
saddle-point state ip s can be continuously shifted (see Eq. (Al)). 

The presence of the zero eigenvalues A^, \f\ and A^ requires some extra efforts in 
evaluating the prefactors in Eq. (22). While the gauge invariance comes from an unphysical 
degree of freedom and hence A^ and Aj 2 ) are simply discarded, the translational invariance, 
however, points to the fact that phase-slip fluctuations in different spatial regions are equally 
probable. The total rate for a transition, say ip n -i ipn, should be obtained by summing 
over all the phase-slip fluctuations across the whole length of the system. This is achieved 
as follows. 

The prefactors £l± in Eq. (22) appear to diverge because of the presence of [A®] -1 / 2 in 
(det Hs)' 1 / 2 . This [A( 3) ]~ 1/2 is contributed by the integral 



2nk B T 



cxp 



A?) 



(25) 



2nk B T V2 

where (l/2)\^[c^] 2 comes from a Taylor expansion at t] s for F[rj], representing the second- 
order term contributed by the component of 77 — t] s in the direction of . Here u^ 3 ) (x) is 



10 



the normalized eigenvector corresponding to and is the coordinate in the direction 
of u( 3) (x): 

H s uf\x) = \^uf\x) : c?> = J dx[uf\x)] T [r,(x) - •,.(*)] ■ 
Given \^ = 0, equation (25) may be rewritten as 



A <V i J^m i d4 , (26) 



2tt V k B T 

It can be shown that / dcf> is an integral proportional to the system length I: J dcf^ = Al, 
and hence 

-i/2 1 laiHl{T) 



A (3) 



2tt V knT 



-Al. (27) 



Thus the prefactors Q± in Eq. (22) are proportional to the system length, as required 
physically by the translational symmetry of the system. In Ref. [3] an analytical expression 
has been derived for A (see Appendix. B). Below we outline a method for evaluating A 
numerically. 

The eigenspace of H s corresponding to the zero eigenvalue is two-dimensional. An or- 
thonormal basis can be constructed from the saddle-point state T) 3 (x). The translational 
invariance of F gives Vt] s {x) as an eigenvector of H s : H s Vrj s (x) = 0. The gauge invariance 
gives u^ 2 )(x) in Eq. (24) as another eigenvector of H s : H s uf\x) = 0. Both Vr) s (x) and 
\if\x) are readily computed from r) s (x) numerically. Based on these two nonorthogonal 
eigenvectors, the eigenvector u^(a;), normalized and orthogonal to uf\x), is obtained: 

u?>(s) = «?> [Vri s (x) - (u?>(s), Vi^s)) uf\x)} , (28) 

where is the normalization factor and (u^(x), Vrj s (x)) stands for the inner product 
/ dx[u^(x)] T Vrj(x). Consider an infinitesimal variation of rj s (x), dr] s (x) = Vr) s (x)dx c , 
where dx c denotes an infinitesimal translation of the phase-slip center. The corresponding 
change of the coordinate in the direction of u^ 3 ) (x) is given by 

dcf ] = (uf\x), dr] s (x)) = (uf\x), Vrj s (x)) dx c = Adx C) (29) 
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where A is related to the normalization factor a( 3 ) in Eq. (28) by A = [a^ 3 )] -1 . Heuristically, 

Eq. (28) defines u^(x) to represent a special direction in the ^-function space. Along 

this direction, the change of r) s (x) is a "pure" translation of the phase-slip center, without 

any rotation of the global phase angle, which is an unphysical degree of freedom. Then in 

Eq. (29), the projection of drj s (x) onto u|, 3 )(x) measures the infinitesimal translation of the 

phase-slip center in the ^-function space, along the physically nontrivial direction of u^(x). 

With the help of Eq. (29) and / dx c = I, Eq. (27) is obtained from Eq. (26). 

To summarize, the dimensionless expressions for the prefactors fl± are obtained as 

|A^| H&AT) [ det/ifn-! 
+ 4tt 2 V k B T [\det"H s \ 
l atH 2 c (T) \ det'H n 

4vr 2 V k B T Af [|det ff ^| 

where the zero eigenvalues and A^ 2 -* are omitted and Eq. (27) is used for [A^ 3 -*] -1 / 2 . Here 
det' in det' H n -i and det' H n indicates that the only zero eigenvalue is to be omitted when 
computing the determinant, and det" in det" H s indicates that the two zero eigenvalues 
are to be omitted when computing the determinant. Using those relevant eigenvectors 
(corresponding to the negative and zero eigenvalues), the matrices H n _i, H n , and H s can 
all be modified into positive definite matrices (as expressed by Eq. (11)), for which the 
determinant ratio can be computed according to Eq. (9). 



1/2 



(30) 



V. NUMERICAL RESULTS 
A. Minimal Energy Path 

The string method has been employed to calculate the MEPs connecting neighboring 
metastable states ip n and ip n -i- All quantities in the numerical calculation are dimensionless. 
The length of the system is / = 327r. The local minima of F are given in Eq. (17) with 
M < 9. (The maximum |n| allowed by \k n \ < k c is 9.) 

We first show in Fig 1 the MEP which connects ip^ to ip^. The string is discretized 
by M = 101 points in the ■0(x)-function space. The initial string is taken from a linear 
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interpolation between tp^x) and ips(x). In order to reach the MEP, the string is evolved 
toward the steady state according to Eq. (3), with the potential force given by 

— — F[rj(x)] = V 2 r; + 77 — rfr). 
or) 

During this process, the string is re-parametrized by arc length every 10 steps. In the 
calculation, ip(x) is represented by a column vector of 2N entries, with the x interval [0,/] 
discretized by a uniform mesh of N = 100 points. Spatial derivatives in the potential force 
are discretized using central finite difference. 

To fix the global rotation of the system, a spring force is applied to the endpoint order 
parameter ^(0). In the form of f = —ICr] 2 (0) with /C = 50, this force restricts ip(0) to the 
real axis. 

The first column in Fig 1 displays a sequence of the configurations along the MEP from 
^4 (x) to ip3(x), and the second column displays the corresponding sequence of Along 
this particular MEP, there is a phase slip of 2ir, nucleated in the middle of the system. 
Through this phase slip, the winding number changes from n = 4 to n - 1=3. From Fig. 
1, it is seen that | ip(x)\ first decreases and reaches zero somewhere (at x = 1/2, see the 
fourth figure from the top), then the phase slip occurs and | '0(^)1 rebounds to accomplish 
the transition. The third figure from the top shows the saddle point ip s (x) between tp^x) 
and i])z{x), which has a locally diminished amplitude and possesses the highest energy along 
the MEP. 

Little [1] first pointed out that a persistent current in a closed loop will not be destroyed, 
"unless a fluctuation occurs which is of such an amplitude that the order parameter is driven 
to zero for some section of the loop". However, the configuration of a vanishing order- 
parameter amplitude somewhere does not necessarily correspond to the lowest saddle point 
between two current-carrying metastable states. Using the stationary Ginzburg-Landau 
equation, Langer and Ambegaokar [2] have obtained the analytical solution for the free- 
energy saddle point. They also pointed out the following: "It is plausible that, from this 
state of locally diminished amplitude, the system will run downhill in free energy through 

13 



a configuration in which the amplitude vanishes somewhere, and finally will achieve the 
configuration in which one less loop in ip occurs across the length L." This picture about 
the transition pathway has been quantitatively confirmed by the MEP obtained here. 

For comparison, we have carried out direct simulations for the motion of ip in the pres- 
ence of thermal noise, using the stochastic equation (20). For k B T <C a£al(T c — T) 2 /(3, 
reasonably clean transition pathways can be obtained from the rare transition events which 
carry the system from one metastable state to the other. Figure 2 displays a sequence of 
ip(x) and IVK^OI) collected along a transition pathway from ij)±{x) to ^(x), calculated for 
I = 327r and ksT = 0.02(x£aQ(T c — T) 2 //3. A comparison of Figs. 1 and 2 shows remarkable 
similarities. The advantage of a MEP is also seen from this comparison: As a smooth path in 
configuration space, the MEP reveals the transition behavior better than those noisy path- 
ways obtained from stochastic simulations. While some fine features of the transition may 
be lost due to the noise in stochastic simulations, they can be well preserved in the MEP. In 
particular, in order to obtain a clean pathway from stochastic simulation, the temperature 
must be kept low enough to reduce local fluctuations, but a low temperature inevitably 
makes the transition events rare and difficult to catch, thus requiring very long simulation 
time. 

Figure 3 shows the energy variation along the MEP from ip^ to ip^. The dimensionless 
free-energy barrier for the transition t/> 4 — ► ip 3 is obtained as AF_ = F[ip s ] — F[ip 4 ] = 0.338. 
Figure 4 shows the energy profile along the MEP from ip 8 to t/>_ 8 . This MEP consists of 16 
segments, each connecting two neighboring metastable states ip n and ip n -i, with n running 
from 8 to -7. 

B. Prefactor 

In calculating the prefactor f2_ in Eq. (30) for n = 4, we use the following procedure: 
(1) From the MEP calculated in Sec. VA, the minimum r\ n and the saddle point r\ s are 
obtained. 
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(2) The (unphysical) degenerate directions at r\ n and u^ 2 ) at rj s in the ^-function 
space are obtained by simple rotation according to Eqs. (23) and (24). 

(3) The degenerate direction uj, 3 ) at r} s in the ^-function space is calculated using Eq. 
(28), and the parameter A defined in Eq. (29) is then obtained to be 0.894. 

(4) The unstable direction at rj s is obtained from the normalized difference be- 
tween two neighboring configurations, evaluated at the saddle point along the MEP. The 
corresponding negative eigenvalue is obtained as = (u^\H s u^) = —0.364; 

(5) The Hessians H n and H s are modified to give two positive definite matrices 



H n = H n + i/W uW * (31) 



T 



and 

3 



H S = H S + Y: u« u« (32) 



1=1 

where v$ and i/W's are all positive parameters; In our calculation the parameters and 
z/W's are set to be 1. 

(6) The ratio det H n / det H s is calculated using Eq. (9). Figure 5 shows the expectation 
value Q(a), defined in Eq. (7) as a function of a in the interval [0, 1]. This interval of a is 
discretized using a non-uniform mesh of 352 points. Since Q(a) varies rapidly near a = 
and 1, more points are distributed near these two ends, with the grid size Aa = 5 x 10~ 4 . 
In the middle of [0,1], a larger grid size Aa = 1.25 x 10~ 2 is used. For each a, the 
stochastic equation (10) is simulated with e = 1, and Q(a) is obtained from a time average 
over 10 5 realizations. The calculated value of the ratio det det if s is 0.608. Using 
this ratio and the computed negative eigenvalue at the saddle point, we obtain the ratio 
det' #4/1 det" H s \ = 1.06. Other sets of values have also been used for and in Eqs. 
(31) and (32), and e in Eq. (7), but the final result of det'if 4 /| det" H s \ is not sensitive to 
those values. 
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C. Rate 



Using I = 32vr and a£o^(T c -T) 2 /[3 = a£H 2 (T) /2tv = 50k B T (used in stochastic simula- 
tion), and the numerical values of AF_ = F[if> s ] - F[ip 4 } = 0.338, = -0.364, A = 0.894, 
and det'-f/zt/l det"if s | = 1.06 obtained in Sees. VA and VB, we are ready to compute the 
rate (in the unit of l/r(T)) 

1 " " V k B T 

for the transition ip4 — > ip^. The exponential factor is e~°- 338x5 ° ps 4.6 x 10~ 8 . The prefactor 
f2_ is evaluated as 

^ = I ~ °^ 64 l x V2n x 50 x 0.894 x 32vr x (1.06) 1/2 w 15.1. 

47T 2 

It follows that the dimensionless rate T_ is approximately 7x 10 -7 . This value is in reasonable 
agreement with what has been estimated through stochastic simulations. 

VI. DISCUSSION 

In this paper we have demonstrated that by using the string method, thermal transi- 
tion rates as formulated in the LAMH theory can be numerically evaluated, even at low 
temperatures. In particular, the pre-exponential factor may also be determined to some 
precision. Thus the "electrical resistance" of a ID superconductor may be evaluated quan- 
titatively. However, it has to be pointed out that quantum tunneling effect, which can be 
important at low temperatures, is not taken into account in the present formulation. Work 
is presently underway to show that quantum tunneling can be similarly treated through 
the string method, thus enabling a complete quantitative account of the current dissipation 
phenomenon in ID superconductors. 
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det' H 4 




\det"H 8 \ 


exp 



27rk B T 



(33) 



APPENDIX A: LANGER-AMBEGAOKAR FREE-ENERGY SADDLE POINT 



By definition, the the lowest saddle point ip s (x) between two neighboring metastable 
states ip n and ip n -i satisfies the stationary Ginzburg-Landau equation (16). Langer and 
Ambegaokar have obtained 



ip s (x) = < \Jl — 3k 2 tanh 



-x 



iV2k s 



(Al) 



where k s is a wave vector determined by the condition 



'1 — 3k 2 

4(1/2) - 0(-Z/2) = kj + 2 tan" 1 ( S 



2irn, 



satisfying k n ^\ < k s < k n . Note that the amplitude of ip s (x) is diminished in a small region 
around x = 0, the phase-slip center. Here we note that translating the phase-slip center 
from x = to x c will produce another saddle point ip s (x — x c ), because the free energy of 
the system is translationally invariant. From the explicit expressions for ip n and ip s , the 
dimensionless energy barriers can be readily obtained: 

iV2 



AF_ 



AF+ = - 
+ 4 



1 

4 

8^2 



, v 1 - 3k 2 s 

l-3k 2 - 8U1 - k 2 ) tan' 1 1—= 



l-3k 2 + 8U1 - ki) TT - tan 



(A2) 



3 v " ^ \ y/2k a 

Here AF_ = F[^ s ]-F[^ n ] and AF+ = F[ip a ] -F[-0 n -i]- Since AF_ < AF+ (for fc s > 0), the 
transition from ■?/>„ to ip n -i is more probable than that from ip n -\ to As a consequence, 
thermally activated phase slips are current-reducing dissipative process. 



APPENDIX B: MCCUMBER-HALPERIN EXPRESSION FOR A 

By writing the real and imaginary parts as the two components of a vector, the saddle- 
point state ip s in Eq. (Al) can be written as 
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cos k s x — sin k s x 

sin k s x cos k s x 

cos k s x — sin 

sin k*x cos 



1 — 3/c 2 tanh 



'1 -3fc s 2 



-V2k s 



(Bl) 



Vr7 s (x) 



-1 






T} s {x) + 


1 





Vfj s (x) 



(B2) 



(B3) 



from which 

cos k s x — sin k s x 
sin /c s x cos /c s a; 

is obtained. This equation indicates that while Vrj s is an eigenvector of H s corresponding to 
the zero eigenvalue, i.e., H s Vf] s = by definition, it can be decomposed into two physically 
distinct components. The first component, 

-1 

1 

is the eigenvector corresponding to the zero eigenvalue \^\ arising from gauge invariance. 
The second component, 

cos k s x — sin k s x 
sin k s x cos k s x 

is the eigenvector corresponding to the zero eigenvalue \f \ arising from translational sym- 
metry. Here a constant A is introduced for normalization. It is easily seen that the global 
rotation of the phase angle is achieved by changing r) s in the direction of while the 
translation of the phase-slip center is achieved by changing rj s in the direction of v^ 3 ). The 

n2 2y/2, 



Vf, s (x)=Av?\x), 



(B4) 



normalization constant A is determined by A 2 = / dx [Vf] s (x)] = 



3A; 2 ) 3 / 2 . The 



change of rj s (x) due to a small change in the location of the phase-slip center, x c , is 



cos k„x — sin k.x 



sin k s x cos k s x 



Vf) s (x)dx c = vi 3 \x)dci?K 



.(3) 



(B5) 



Substituting Eq. (B4) into Eq. (B5) yields dcf^ = Adx c and thus / dc^ = A / dx c = Al. 
Therefore the sample length dependence arises naturally from the translational degeneracy. 
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We want to point out that in general, vf\x) is not orthogonal to u^(x). From the 
inner products 











T 








-1 






-1 


Jdx 




1 


r) s (x) 




1 



T] S {X) 



j dx ir) s (x)] T Vs (x) = (1 - kl)l - 2^1 - 3kl 



(B6) 



/ dx 



cos k q x — sin k.x 



sin k s x cos k s x 





T 












cos k s x 
sin k s x 


— sin k s x 




Vfj s (x) 




cos /c s x 


Vf] s (x) 



2V2 



Jdx[Vf, s (x)] T Vfi s (x) = J^-(l-3k 2 s f 2 , 



(B7) 









T r 




-1 


r} s (x) 






1 







2v^fc sV /l-3fc2, 



cos keX — sin /c,a; 



sin k s x cos /c s x 



Vfj a (x) 



we obtain 



/ 



dx 



uf\x) 



T 



vf\x) 



1/2 



(B8) 



(B9) 



(1 - fc2)Z _ 2^2^ 1 - 3A;2j (1 - 3k*)V* 
which approaches zero as k s — > and/or / — > 00. So the two are orthogonal only in these 



limits. 
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FIGURES 




FIG. 1. Minimal energy path from ip^ to "03- Left column: (Reip, Imip) as a function of x; 
Right column: \ip\ as a function of x. The figures at the top and bottom correspond to ip4 and ^3) 
respectively. The third from the top is for the saddle point. In the fourth figure, reaches zero 
in the middle. 
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FIG. 2. Transition pathway obtained from stochastic simulation, from tp^ to ip^ with thermal 
noise. Left column: (Reip, Imip) as a function of x; Right column: as a function of x. The 
figures at the top and bottom correspond to ^4 and V>3> respectively. The third from the top is 
closest to the saddle point. In the fourth figure, |^| reaches zero in the middle. 



22 



-21.5 




23 



CO 



-12 
-14 
-16 
.-18 
-20 
-22 
-24 
-26 



n=-8 


I 


1 

n=8 


-7 




7 


-6 | 




6 




-5 5 






-4 4 






-3 \A . , |VJ 3 
-1 Q 1 


1 



-1 







1 



FIG. 4. Scaled free energy F evaluated along the MEPs connecting a sequence of metastable 
(n 7^ 0) and stable (n = 0) states, plotted as a function of the arc length s in the ^(x)-function 
space. The scaled system length is I = 32tt and the zero-current state (n = 0) is taken as the 
reference point where F = Sir and s = 0. The maximum |n| allowed by \k n \ < k c is 9. The plot 
here is up to n = ±8. 
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